In [ ]:
import io
import math
import os
import pathlib
import holoviews as hv
import hvplot.pandas
import numpy as np
import pandas as pd
import pymap3d
import xarray as xr
hv.extension("bokeh")
np.set_printoptions(suppress=True)
In [ ]:
def get_skews_and_base_cfls(lons, lats, depths) -> np.ndarray:
# The shape of each one of the input arrays needs to be (3, <no_triangles>)
#ell = pymap3d.Ellipsoid.from_name("wgs84")
ell = pymap3d.Ellipsoid(6378206.4, 6378206.4, "schism", "schism")
local_x, local_y, _ = pymap3d.geodetic2enu(lats, lons, depths, lats[0], lons[0], depths[0], ell=ell)
areas = (local_x[1] * local_y[2] - local_x[2] * local_y[1]) * 0.5
rhos = np.sqrt(areas / np.pi)
max_sides = np.maximum(
np.sqrt(local_x[1] ** 2 + local_y[1] ** 2),
np.sqrt(local_x[2] ** 2 + local_y[2] ** 2),
np.sqrt((local_x[2] - local_x[1]) ** 2 + (local_y[2] - local_y[1]) ** 2),
)
skews = max_sides / rhos
base_cfls = np.sqrt(9.81 * np.maximum(0.1, depths.mean(axis=0))) / rhos / 2
return skews, base_cfls
def get_skews_and_base_cfls_from_path(path: os.PathLike[str] | str) -> np.ndarray:
ds = xr.open_dataset(path, engine='selafin')
tri = ds.attrs['ikle2'] - 1
lons = ds.x.values[tri].T
lats = ds.y.values[tri].T
depths = - ds.B.isel(time=0).values[tri].T
skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)
return skews, base_cfls
In [ ]:
file = "/home/tomsail/Documents/work/models/meshes/slf/v1p2.slf"
file = "/home/tomsail/Documents/work/models/meshes/slf/v2p1.slf"
ds = xr.open_dataset(file, engine='selafin')
skews, base_cfls = get_skews_and_base_cfls_from_path(file)
In [ ]:
CFL_THRESHOLD = 0.4
for dt in (1, 50, 75, 100, 120, 150, 200, 300, 400, 600, 900, 1200, 1800, 3600):
violations = (base_cfls * dt < CFL_THRESHOLD).sum()
print(f"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%")
1 5887124 100.00% 50 1566579 26.61% 75 1095580 18.61% 100 824343 14.00% 120 664172 11.28% 150 344840 5.86% 200 158392 2.69% 300 61959 1.05% 400 11172 0.19% 600 3120 0.05% 900 1934 0.03% 1200 522 0.01% 1800 11 0.00% 3600 0 0.00%
In [ ]:
pd.DataFrame({"skew": skews}).describe()
Out[Â ]:
| skew | |
|---|---|
| count | 5.887124e+06 |
| mean | 2.826159e+00 |
| std | 1.493293e-01 |
| min | 2.506716e+00 |
| 25% | 2.724916e+00 |
| 50% | 2.790498e+00 |
| 75% | 2.894929e+00 |
| max | 2.021513e+01 |
In [ ]:
df = pd.DataFrame({"cfl": base_cfls * 400})
df[df.cfl < 0.4].describe()
df[df.cfl < 0.4].hvplot.hist()
Out[Â ]:
In [ ]:
tri = ds.attrs['ikle2'] - 1
nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=["lon", "lat", "depth"])
elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=["no_nodes", "n1", "n2", "n3"])
elements = elements.assign(base_cfl=base_cfls)
elements.head()
Out[Â ]:
| no_nodes | n1 | n2 | n3 | base_cfl | |
|---|---|---|---|---|---|
| 0 | 3.0 | 802081.0 | 65536.0 | 2267476.0 | 0.002489 |
| 1 | 3.0 | 585525.0 | 65536.0 | 802081.0 | 0.002548 |
| 2 | 3.0 | 65536.0 | 585525.0 | 780858.0 | 0.002697 |
| 3 | 3.0 | 251963.0 | 65536.0 | 780858.0 | 0.002613 |
| 4 | 3.0 | 65536.0 | 251963.0 | 820904.0 | 0.002469 |
In [ ]:
min_cfl_per_node = pd.concat([
elements[["n1", "base_cfl"]].groupby(["n1"]).base_cfl.min(),
elements[["n2", "base_cfl"]].groupby(["n2"]).base_cfl.min(),
elements[["n3", "base_cfl"]].groupby(["n3"]).base_cfl.min(),
], axis=1).min(axis=1)
min_cfl_per_node.head()
Out[Â ]:
1.0 0.005506 2.0 0.007554 3.0 0.003865 4.0 0.018533 5.0 0.021136 dtype: float32
In [ ]:
dt = 200
df = nodes.assign(
cfl=min_cfl_per_node * dt,
# CFL_violation nodes have a value of 1 if there is no violation and 4 if there is a violation.
# We do this in order to plot the points with a different size
cfl_violation=((min_cfl_per_node * dt < CFL_THRESHOLD) * 3) + 1
)
df.head()
Out[Â ]:
| lon | lat | depth | cfl | cfl_violation | |
|---|---|---|---|---|---|
| 0 | -180.000000 | 90.000000 | -4228.0 | NaN | NaN |
| 1 | 134.560318 | -65.732597 | -107.0 | 1.101187 | 1.0 |
| 2 | -13.422298 | 29.300503 | -68.0 | 1.510851 | 1.0 |
| 3 | 12.516141 | 54.809700 | -20.0 | 0.772975 | 1.0 |
| 4 | 131.169098 | -0.395918 | -300.0 | 3.706690 | 1.0 |
In [ ]:
plot = df[df.cfl_violation == 4].hvplot.points(
'lon',
'lat',
c="depth",
cmap="jet",
geo=True,
tiles="EsriImagery",
).options(
width=1200, height=900
)
len(df[df.cfl_violation == 4])
Out[Â ]:
135053
In [ ]:
plot
Out[Â ]: